module sortff
    implicit none
    private

    public :: sort, sorted, sorted_order

    interface sorted_order
        module procedure sorted_order_int
        module procedure sorted_order_sp
        module procedure sorted_order_dp
    end interface

    interface sorted
        module procedure sorted_int
        module procedure sorted_sp
        module procedure sorted_dp
    end interface

    interface sort
        module procedure sort_int
        module procedure sort_sp
        module procedure sort_dp
    end interface
contains
    pure recursive function sorted_order_int(array) result(sorted_indices)
        integer, intent(in) :: array(:)
        integer, allocatable :: sorted_indices(:)

        logical, allocatable :: greater_than_pivot(:)
        integer :: i
        integer, allocatable :: indices(:)
        integer, allocatable :: indices_equal_pivot(:)
        integer, allocatable :: indices_greater_than_pivot(:)
        integer, allocatable :: indices_less_than_pivot(:)
        logical, allocatable :: less_than_pivot(:)
        integer :: n
        integer :: pivot
        integer, allocatable :: sorted_greater_than(:)
        integer, allocatable :: sorted_less_than(:)

        n = size(array)
        select case (n)
        case (0)
            allocate(sorted_indices(0))
        case (1)
            sorted_indices = [1]
        case (2)
            if (array(1) > array(2)) then
                sorted_indices = [2, 1]
            else
                sorted_indices = [1, 2]
            end if
        case default
            pivot = (n/2 + 1)
            allocate(indices, source = [(i, i = 1, n)])
            less_than_pivot = array < array(pivot)
            greater_than_pivot = array > array(pivot)
            indices_less_than_pivot = pack(indices, less_than_pivot)
            indices_greater_than_pivot = pack(indices, greater_than_pivot)
            indices_equal_pivot = pack(indices, .not.(less_than_pivot.or.greater_than_pivot))
            allocate(sorted_less_than, source = sorted_order(array(indices_less_than_pivot)))
            allocate(sorted_greater_than, source = sorted_order(array(indices_greater_than_pivot)))
            sorted_indices = &
                    [ indices_less_than_pivot(sorted_less_than) &
                    , indices_equal_pivot &
                    , indices_greater_than_pivot(sorted_greater_than)]
        end select
    end function

    pure recursive function sorted_order_sp(array) result(sorted_indices)
        real, intent(in) :: array(:)
        integer, allocatable :: sorted_indices(:)

        logical, allocatable :: greater_than_pivot(:)
        integer :: i
        integer, allocatable :: indices(:)
        integer, allocatable :: indices_equal_pivot(:)
        integer, allocatable :: indices_greater_than_pivot(:)
        integer, allocatable :: indices_less_than_pivot(:)
        logical, allocatable :: less_than_pivot(:)
        integer :: n
        integer :: pivot
        integer, allocatable :: sorted_greater_than(:)
        integer, allocatable :: sorted_less_than(:)

        n = size(array)
        select case (n)
        case (0)
            allocate(sorted_indices(0))
        case (1)
            sorted_indices = [1]
        case (2)
            if (array(1) > array(2)) then
                sorted_indices = [2, 1]
            else
                sorted_indices = [1, 2]
            end if
        case default
            pivot = (n/2 + 1)
            allocate(indices, source = [(i, i = 1, n)])
            less_than_pivot = array < array(pivot)
            greater_than_pivot = array > array(pivot)
            indices_less_than_pivot = pack(indices, less_than_pivot)
            indices_greater_than_pivot = pack(indices, greater_than_pivot)
            indices_equal_pivot = pack(indices, .not.(less_than_pivot.or.greater_than_pivot))
            allocate(sorted_less_than, source = sorted_order(array(indices_less_than_pivot)))
            allocate(sorted_greater_than, source = sorted_order(array(indices_greater_than_pivot)))
            sorted_indices = &
                    [ indices_less_than_pivot(sorted_less_than) &
                    , indices_equal_pivot &
                    , indices_greater_than_pivot(sorted_greater_than)]
        end select
    end function

    pure recursive function sorted_order_dp(array) result(sorted_indices)
        double precision, intent(in) :: array(:)
        integer, allocatable :: sorted_indices(:)

        logical, allocatable :: greater_than_pivot(:)
        integer :: i
        integer, allocatable :: indices(:)
        integer, allocatable :: indices_equal_pivot(:)
        integer, allocatable :: indices_greater_than_pivot(:)
        integer, allocatable :: indices_less_than_pivot(:)
        logical, allocatable :: less_than_pivot(:)
        integer :: n
        integer :: pivot
        integer, allocatable :: sorted_greater_than(:)
        integer, allocatable :: sorted_less_than(:)

        n = size(array)
        select case (n)
        case (0)
            allocate(sorted_indices(0))
        case (1)
            sorted_indices = [1]
        case (2)
            if (array(1) > array(2)) then
                sorted_indices = [2, 1]
            else
                sorted_indices = [1, 2]
            end if
        case default
            pivot = (n/2 + 1)
            allocate(indices, source = [(i, i = 1, n)])
            less_than_pivot = array < array(pivot)
            greater_than_pivot = array > array(pivot)
            indices_less_than_pivot = pack(indices, less_than_pivot)
            indices_greater_than_pivot = pack(indices, greater_than_pivot)
            indices_equal_pivot = pack(indices, .not.(less_than_pivot.or.greater_than_pivot))
            allocate(sorted_less_than, source = sorted_order(array(indices_less_than_pivot)))
            allocate(sorted_greater_than, source = sorted_order(array(indices_greater_than_pivot)))
            sorted_indices = &
                    [ indices_less_than_pivot(sorted_less_than) &
                    , indices_equal_pivot &
                    , indices_greater_than_pivot(sorted_greater_than)]
        end select
    end function

    pure function sorted_int(array) result(sorted)
        integer, intent(in) :: array(:)
        integer, allocatable :: sorted(:)

        sorted = array(sorted_order(array))
    end function

    pure function sorted_sp(array) result(sorted)
        real, intent(in) :: array(:)
        real, allocatable :: sorted(:)

        sorted = array(sorted_order(array))
    end function

    pure function sorted_dp(array) result(sorted)
        double precision, intent(in) :: array(:)
        double precision, allocatable :: sorted(:)

        sorted = array(sorted_order(array))
    end function

    pure subroutine sort_int(array)
        integer, intent(inout) :: array(:)

        array = sorted(array)
    end subroutine

    pure subroutine sort_sp(array)
        real, intent(inout) :: array(:)

        array = sorted(array)
    end subroutine

    pure subroutine sort_dp(array)
        double precision, intent(inout) :: array(:)

        array = sorted(array)
    end subroutine
end module
